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Abstract 

There is a growing interest in developing covariance functions for processes on the surface of a 
sphere due to wide availability of data on the globe. Utilizing the one-to-one mapping between 
the Euclidean distance and the great circle distance, isotropic and positive definite functions in 
a Euclidean space can be used as covariance functions on the surface of a sphere. This approach, 
however, may result in physically unrealistic distortion on the sphere especially for large distances. 
We consider several classes of parametric covariance functions on the surface of a sphere, defined 
with either the great circle distance or the Euclidean distance, and investigate their impact upon 
spatial prediction. We fit several isotropic covariance models to simulated data as well as real 
data from NCEP/NCAR reanalysis on the sphere. We demonstrate that covariance functions 
originally dehned with the Euclidean distance may not be adequate for some global data. 

Keywords: Covariance function; Euclidean distance; Geopotential height; Great circle 
distance; Process on a sphere. 


1. Introduction 

In geophysical and environmental sciences, data often come in a global scale. Covariance 
functions for global data sets need to be positive definite on the surface of a sphere and the 
distance computation is important in spatial modeling. For an integer d > 1, dehne = {x G 
]^d+i . ii^ii _ ^ (d-dimensional) sphere with radius r, where ||x|| is the Euclidean norm 

of X G We also define the great circle distance on S!^ by 6*(x, y) = r x arccos((x, y)) where 

(•, •) denotes the usual inner product on 

We consider the surface of the Earth as the spatial domain. Let denote the surface of 
the Earth, where R denotes the Earth’s radius (we approximate the Earth as a perfect sphere). 
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The chordal distance between the two points, Si = (Li,/i) and S 2 ={L 2 ,l 2 ), on iS|. (L and / denote 
latitude and longitude, respectively) is given by 

ch(si, S 2 ) = 2i?[sin^{(Li — L 2 )/ 2 } + cosLicosL 2 sin^{(/i — / 2 )/ 2 }]^^^. 


The great circle distance between the two locations then is given by 6* = 0 ( 81 , 82 ) = 2R x 
arcsin{ch(8i, 82 ) 7 ( 2 /?)}. The chordal distance is simply the Euclidean distance penetrating the 
spatial domain on th e surface of the Earth, and producing a straight line approximation to the 


g reat circle c 


istance flBaneried . 


200 ^. 


Yadrenkol (119831 ) pointed out that any covariance function in can be considered as a covari¬ 


ance function for processes on using the chordal (i.e . Euclidean) dis tance. This con struction 
can p rovide a rich class of covariance functions on (IGneitind. Il9991) . As argued in iGneiting 
(120131 ). the great circle distance is a physically most natural distance metric for processes on a 
sphere. However, literature on covariance modeling using the great circle distance on the surface 
of a sphere is scarce due to its mathematical challenge. Some efforts have been made in examining 


the validity of severa l parametric covariance functions on the surface of a sphere (iHnang et ah 


2011 


Gneitingl . l2013l ) and in developing valid para metric covariance functions with the 


cle distance from various construct i onal approaches 


2013 


Guinness and Fuentesl. 


Although 


Huang et ah 


2013; 


Jeong and 


(120 llh and 


Gneiting 


un 


(Du and Ma. 

2012; 

Du et ah. 

2013: 


great cir- 


Gneiting . 


20151) . 


(120131) studied validity of covariance functions 


dehned with either the great circle distance or the Euclidean distance in details, th e impact 


upon 


parameter estimation and prediction has not been studied well. According to 


Baneriee 


(120051) . careless formulation of distances can lead to poor prediction with wrong estimation 


of the spatial range. Note that this study considered the Matern covariance model using the 
great circle distance, which may not be positive dehnite on the s urface of a sph ere, unless the 


2001 


Gneitingl. 120131) . 


smoothness parameter, z/, is z/ G (0,0.5] (IMiller and Samko . 

In this paper, we consider several positive dehnite functions on and for d = 1,2, and 
compare them in simulation studies and real data applications.The rest of the paper is organized 
as follows. In Section 121 we discuss some characteristics of covariance functions on a sphere. 
Then we present two simulation studies on and in Section l3l Section 0] illustrates real 
application results to geopotential height data set. Finally, Section 0] concludes the paper with 
discussion. 
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2. Characteristics of covariance functions on a sphere 


We first review some known results on covariance functions in the Euclidean space as well as 
those on the surface of a sphere. A function / : x —)■ M is called positive definite if 


E"i=iCiCj/(xi,Xj) > 0 


( 1 ) 


for all hnite n G N, all distinct points xi,... ,x„ G and all real Ci,... ,c„. A function is 
strictly positive definite when the inequality in ([T]) is strict (unless ci = C 2 = • • • = = 0). 

For a real random field Z in with E{Z{x.)}‘^ < oo for all x G the covariance function 
A(x, y) = cov{Z(x), Z{y)} must satisfy the condition in (1). The random field Z is called weakly 
stationary if its mean function is constant, it has finite second moments, and its covariance 
function can be written as cov{Z(x), Z(y)} = A'(x — y) for all x, y G M'^, and a positive 
definite function K, i.e. the covariance function of Z, depends on x and y only through x — y. 
Furthermore, if its covariance function satishes cov{Z(x), Z(y)} = (p(||x — y||) for a positive 
definite function (p, then the random held Z is weakly isotropic. An isotropic prope rty for 


proce sses in can be thought as an invariance property under translation and rotation fISteinl . 


1999h 


The covariance function of a random held and the smoothness of its realization are related to 
mean square properties of the random held. The random held Z is called mean square continuous 
at X if 

E{Z{y) - Z(x)}^ -^0 as y -> X. 

For weakly stationary random held, mean square continuity is equivalent to the fact that the 
covariance fu nction is continuous at the origin, but it does not imply continuity of its realization 
(jSteinl . 1199911 . Moreover, a random held Z on M with hnite second moments is mean square 
diherentiable at t if there exists Z'{t) = lim„^oo{Z(t + h„) — Z{t)}/hn in for sequences 
hness of a rand om held can be determined through the number of mean 


0. The smoo 

square derivatives. iGneitind f)2013l) dehned the class of with the correlation functions of 
mean square continuous, stationary and isotropic random helds in Every positive dehnite 
function (p : [0, oo) —)■ M with </?(0) = 1 is the correlation of an isotropic process and the members 


of $2 anc 


(IGneiting 


$3 ar e the cornerstones for covariance models for spatial data in a planar domain 


20131) . An isotropic property on a sphere means the covariance function depends 


on distance only. That is, a random held Z on 5^ is called isotropic if its covariance function 
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satisfies cov{Z(x), Z(y)} = '0(6'(x, y)) for all x, y G S^. We then similarly define \l/d, the class 
of continuous, isotropic covariance functions : [0, tt x r] ^ M on 

Since a sphere can be viewed as a subset of the Euclidean space, valid covariance functions on 


X can be restricted to x when the Euc 
chordal distance on 


idean distance is used (equivalently, the 


Yadrenkd fll983l) and lYagloml fll987l) pointed out that if (yC is a member 


of the class <h(i+i, then the function y9[2r sin{0/(2r)}], with the Euclidean distance expressed 
in terms of great circle distance as 2r sin{0(x,y)/(2r)}, belongs to the class tl'd- Since there 


are vario us po s itive 


families (ISteinl. 


1999 


dehnite function s, including the Matern class and the generalized Cauchy 


Gneitind. 120131 ) that are isotropic covariance functions for processes in 


this mapping from (p G $3 to "0 ^ provides a useful way to construct a rich parametric class 
of isotropic covariance functions on S^. This mappi ng preserv e s the interpretation of parameters 


such as scale, range, smoothness, and fractal ind ex flGneitind. 


It 


Guinness and Fnentes 

f2Q13 

), 

Jeong and Jnn 


(120151 )) that when Matern class with the Euclidean distance and that with the great circle dis¬ 


tance are compared in terms of model £t and prediction, often Matern model with the Euclidean 
distance performs better. This may be due to th e restriction on the sm oothness parameter for 
the Matern class with the great circle distance. IJeong and JuiJ (120151) proposed a method to 
overcome such limitation on the smoothness parameter for the Matern class with the great circle 
distance, but they found that the Matern class with the Euclidean distance is equivalent or often 
better compared to the covariance models specifically developed for processes on the sphere. Our 
goal in this paper is to study cases that are not previously considered in the literature, and to 
explore cases where there are significant differences (improvements) of covariance models defined 
on the sphere as opposed to the covariance models projected from the Euclidean space. 

We focus on the fact that there are some fundamental differences between covariance models 
originally dehned on the surface of a sphere and those in the Euclidean space. For instance, there 
exists a lower bound on isotropic correlation function in the Euclidean space. A function (p is an 
isotropic correlation function in if and only if it has the form, ip{t) = Aii{tu)dG{u), where 
dG{u) = 1 and G is nondecreasing. Note that Ad(r) = J(rf_ 2 )/ 2 (?") 

where J is a Bessel function. Then, for all t, 


(fit) > inG>oAd(s) 


( 2 ) 


(jSteinl . I 1999 I) . This implies that valid correlation functions on constructed through the map- 
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ping described above from (p G $3 cannot have values less than infs>os ^sins = —0.218. In 
particular, the Matern class yields nonnegative correlations only. Although the importance of 


the Matern family is highlighted by 


SteinI f)l999l ) because of its flexibility with regard to the local 


behavior of the processes, it might not be appropriate in applications where there is signihcantly 
negative spatial correlation. In fact, many of the isotropic covariance functions in used in the 
literature take non-negative values only. 

We also focus on the fact that on the sphere, correlation between two points large distance apart 
may not necessarily be small (compared to correlation between nearby two points). In fact, if 
there is a wave traveling around the sphere, two points nearly maximum possible distance apart 
may be perfectly positively (or negatively) correlated, which cannot happen in the Euclidean 
space. 

We now list several parametric classes of covariance functions defined on the surface of a sphere, 
or dehned originally in the Euclidean space (then can be used through the projection). Some of 
the models are used in the simulation and data examples. We only consider isotropic covariance 
models on a sphere. Functions in the class d'rf are characterized in terms of an infinite sum 


of Gegenbauer po! 


(Schoenberg. 

1942; 

Gneiting. 

2013) 


ynomials with n onnegative coefficients and cosine of the great circle distance 


Gneitind. l2013h . For d > 1, the class consists of the functions of the form 


m = « € [O.TT X r], 


degree n, Cn 


(Schoenberg. 

T). d oUl 

1942; 

.11 / 

Ghen et ah. 

fl ^n,d ■ 

2003) 


bn d = 1 and the Gegenbauer polynomial of 


20031) . Moreover, the class \hoo consists of the 


functions with the following form 

= E“=o^n(cos(d/r))”, d e [0,7r X r]. 


(3) 


with nonnegative coefficients bn such that = 1- The infinite sum is strictly positive 

dehnite on Sf when the coefficients bn,d and bn are strictly positive for inhnitely many odd and 
infinitely many even integers n, and only a few closed forms, such as the multiquadra tic family, 
for su ch infinite sums are known in general. The multiquadratic covariance function (jGneiting . 
2013 1 is defined by 


= (T^(l — /{I -f — 2r cos(d/r)}'^, 6 * G [0, vr x r] 

from a standard Taylor series of (jS]), when > 0, c > 0, and r G (0, 1). The Matern class, given 


5 








































as 


ip{t) = ^r(z/) ^{t/aYlCy{t/a), t>0, 


( 4 ) 


where the parameters, a‘^,a,i> > 0, are marginal variance, spatial range, and smoothness param¬ 
eters , respectively, is positive dehnite in for any d G N with the Euclidean distance fjSteinl . 
1999ri. Here, JC,, i s the modihed Bessel function of the second kind of order u. 


Gneitingi (120131 ) showed that completely monotone functions (that have derivatives of all 


orders with (—> 0 for all nonnegative integers k and all positive t) including the 
power exponential, Matern, and generalized Cauchy families are positive dehnite, through the 
restriction of a function ip : [0, cxo) —)■ M to the interval [0 , tt x r]: ip = 9 ^[o, 7 rxr] under applicable 
conditions, on S!^ of any dimension. One necessary condition of the membership in the class 
'hrf is that either the fractal index or the smoothness parameter requires to satisfy (3 G (0, 1] or 
z/ G (0,0.5], respectively. The powered exponential family dehned by 

YiO) = (J^exp[—{ 6 */(cr)}4], d G [0, tt x r] 


where > 0 , c > 0 i s valid on any c 


z/ G (0, 0.5] similarly. I.Teong and ,Tnn 


i mensi onal ii (3 E (0,1], and the Matern family requires 
(120151 ) compare Matern class with great circle distance to 


that the Euclidean distance for spatial data on the surface of a sphere. 

Compactly su pported members of the class $3 may be valid on through Y = </ 5 [o, 7 rxr 


(jCneitind. l2013l) . The spherical and the Wendland’s functions on $ 3 , remain valid with direct 


substitution of the Euclidean distance b y the great circle distance on d = 1, 2, 3. In particular, 
the C^-Wendland covariance function (IWendlandl . Il995l) . dehned as 


^{0) = a {1 + {dT)/{cr) + 0 {t - l)/(3c r )}{1 - d/(cr)}+, d G [0,7r x r]. 


( 5 ) 


where c G (0, tt] is a support parameter and r > 6 is a shape parameter, has 4 derivatives at the 
origin and thus may be suitable for smooth data on 
The following covariance function (dehned in the Euclidean space) models a hole ehect. 


ip{t) = a‘^{a/t) sm(t/a), t > 0, 


( 6 ) 


with 9 ?( 0 ) = (T^, and it is called the wave covariance function. This function may be deal with the 
situation where correlation between two points far apart may have bigger (in magnitude) than 
correlation between two points closer, or oscillating pattern in the correlation fun ction. This 
belongs to the class for d = 1 , 2 , 3, but it is not valid on Sp (IHuang et ahl. I 2 OIII) . Thus, this 
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function can only be applied to the data on the sphere through the projection from the Euclidean 
s pace. _ 


Schoenberg! (119421 ) noted that the class 4/^ enjoys the useful closure properties. The class 


is convex, closed under products, and closed under limits which are continuous. For examples, if 
e Td, then + (1 - G for every A G [0,1] and x ip 2 {d) G 

Moreover, if 'ipniO) G 'ipniO) —> '^{9) as n — )■ oo, and ■^(6') is continuous, then also '0(^) £ "^d- 
These closure properties of the class offer additional flexibility to model negative correlations. 
If we use convex sums or products of parametric families of correlation functions on S!^ in terms of 
the great circle distance with a Legendre function of the form i/j{ 9) = cos(9/r) 


inclu ding the case cos(6'/r) for d = 1,2,3, we can easily model negative correlations flGneitine . 
291311 . 

The sine-power covariance function is dehned by 


2^(9) = - (sin 0 G [0, tt x r]. 


(7) 


where the parameter f3 G (0,2] corresponds to the fractal index and controls the smoothness of 
the process. This covariance function operates directly on a sphere (belongs to the class for 
all dimensions). On the other hand, the cosine function, ^p{t) = cos{nt), where t > 0, belongs to 
the class <hi only for any value of n > 0. On S!^, 'ip{9) = cos{n9/r) for 0 G [0, tt x r], is non-strictly 
positive dehnite for all dimensions when n = 1. For non-integer values of n > 0, 'ip{9) is not valid 
on S} and fo r integ er values of n > 2, it is non-strictly positive dehnite only on but not on 
( Gneitine . 2013 1. 


When there is no signihcant negative correlations, and when the spatial range is small, we do 
not expect that using the great circle distance or the Euclidean distance may make signihcant 
diherence. Moreover, when the prediction location is surrounded by enough number of estimation 
locations, we do not expect signihcant diherence between models dehned on the sphere and 
models dehned in the Euclidean space (and used through projection). Therefore, for the rest of 
the paper, we mainly consider the situations where either spatial range is large, or when there are 
not many number of estimation locations close to the prediction locations. Such situation may 
arise in real application, for example, when we are interested in predicting over the ocean with 
most observations on the land, or when we are interested in predicting near the poles with not 
many observations near the poles. We also consider cases where there is signihcantly negative 
correlations at large distances. 
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3. Simulation studies 


We present two simulation studies on Sj and S^. In the first example (Section 13.11) . the 
truth is generated from an exponential covariance model with the great circle distance. In the 
second example fSection 13.21) . th e truth is generated from oscillating Matern covariance model 
implemented in R-INLA package flMartins et ahl . 120131 ). 


3.1. Example on Si 

We simulated mean zero Gaussian random helds on Si using exponential covariance mod¬ 
els with various range parameters using the great circle distance. The exponential covariance 
function on S^! is dehned by 

■^(6*) = a^exp(—3/a), 9 E [0,7r x r], 


where > 0 is a variance and a > 0 is a spatial range parameter. We then compared htted 
results using the true covariance model as well as exponential covariance model with the Euclidean 
distance, i.e. (p(t) = (T^exp(—t/a), t > 0. Exponential covariance functions belong to both $2 
and Ti, so they are valid on the surface of a sphere regardless of distance. We set the marginal 
variance = 1 and considered the various spatial ranges a = 27r, I.Stt, tt, vr/l.S, 7r/2,7r/4. We 
randomly selected 100 locations (angles) for parameter estimation from (7r/2,37r/2) on a unit 
circle, and 10 hxed and equally spaced locations for prediction from [0,7r/2) as in Figure [H We 
compared the covariance models in terms of prediction using the root m ean square error (RMSE) 
as well as the two popular scoring rules (IGneiting and Raftervl. 120071) . the mean absolute error 
(MAE) and the continuous ranked probability score (GRPS). We repeated this experiment 100 
times: sampling locations are different each time and prediction locations are hxed. 

Figures [2](a)-(d) present MAE and GRPS from the two models, the exponential models using 
great circle and chordal distances, displayed against 10 hxed prediction locations, for the true 
spatial range values are 2n and 7r/4. For the larger spatial range, a = 2 tt, we observe that 
the exponential model using the great circle distance performed signihcantly better than that 
using the chordal distance. Except for prediction locations that are relatively close to sampling 
locations, there are considerable diherences between two models in prediction errors. On the 
other hand, for the smaller spatial range, a = 7r/4, there is no signihcant diherence b e tween two 


models in both predic tion errors, which agrees with hndings in 


Guinness and FuentesI (120131) and 


Jeong and JunI (120151) . 


















3.2. Example on S]^ 

We considered mean zero Gaussian random fields on the surface of the Earth (with radius 
R = 6, 371 (km)) with a Matern covariance function from oscillating stochastic partial differential 
equations (SPDE) models with the spdel class from the R-INLA version 0.0-1413638221. We set 
(T = 1, K = 0.5, T = KX V = 2 and Oosc = 0.3. Here k is the spatial scale parameter, 

r controls the variance with a, v controls the smoothness of the process, and Oosc controls the 
strength of oscillation. We subtracted the constant mean (average over all locations) to have 
mean zero residual fields as in Figure [3](a). There are 128 longitude points and 64 latitude points, 
and the size of the data is 8,192. We randomly selected 300 locations where values are smaller 
than 0 for parameter estimation, and 100 locations where values are larger than 1 for prediction. 
We repeated this procedure 100 times: all locations are different each time. It is clear from 
Figures [ni(a) and (b) that values at large distance lags are negatively correlated. We compared 
fitted results from a Matern covariance model with the Euclidean distance (MC) to those from 
a Matern covariance model with the great circle distance (MG). Moreover, to deal with negative 
correlations, we use convex sum of valid covariance functions with the great circle distance (G): 

iIj{6) = o-^[A{l — (sin ^)^} -|- (1 — A) cos(6'/i?)], jS e (0, 2] and A G (0,1), 

We also considered the hole-effect model with the Euclidean distance (H) defined in (E]), a model 
defined with the Euclidean distance, for a comparison. 

Table [T] contains results of parameter estimation and prediction for the various models consid¬ 
ered. From Table [H we observe that both MG and MG have large estimates of the spatial range 
parameter and MG has smaller prediction errors than MG. Although MG has better fit than MG 
in terms of maximum log-likelihood values, MG leads to poor prediction, possibly due to the large 
estimate of the spatial range. Note that the best model in terms of prediction errors is G. For 
G, the estimate of A is close to 1 and resulting estimated model is dominated by the sine-power 
model. Nevertheless, its correlation function allows much smaller correlation values for large dis¬ 
tance lag, compared to the models for the Euclidean space. On the other hand, for H, although 
it allows negative correlations unlike Matern covariance model, it resulted in poor model fit and 
spatial prediction. From Figures ID) a)-(d), we observe significant differences of prediction errors 
between models using the great circle and Euclidean distances for prediction locations who are 
relatively far away from their nearest sampling locations. Overall, G outperforms MG and MG 
in prediction. 
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4. Application 


4 . 1 . Data and mean structure 

We consider geopotential height data on a global scale. The geopotential height approximates 
the actual height of a surface pressure at certain level above mean sea level. It is an adjustment 
to geometric height using the variation of gravity with elevation and latitude. The study of the 


geo potential height migh ; 


the 


3 e im portant in learning abnormal weather phenomena. According to 


Hafez and Almazronil (120141 ). the geopotential height at level 500 hPa plays a dominant role 


in controlling weather and climate conditions. Moreover, it became evident that the variability 


of global geopotentia 


last several decades (IMarshalll . 


height is clear l y impacted by g. 


2002 


Zhu et ah 


2002 


obal warming and climatic indices over th e 


Hafezl. 


2012 : 


Hafez and Almazronil. 120141) . 


Data sets are obtained at level 500 hPa from NCEP/NCAR reanalysis project (http: //www. esrl .noaa. gov 
The NCEP/NCAR reanalysis 1 project uses a state of the art analysis /forecast syst e m to per¬ 


form data assimilation using past data from 1948 to the present (see iKalnav et ahl fjl996l) for 


more detailed information). The output values are given on regular grids, and there are 144 
longitude points and 73 latitude points (the size of the data is 10,512). We used Boreal summer 
geopotential height in the northern hemisphere (June, July, and August; JJA), and for each grid 
we computed a pointwise mean as the average over 2014. The unit is meter for the geopotential 
height and kilometer for distances. For variance stabilization, we took a square root transform 
of the data. 

We decompose the data into its mean structure (large scale variation) and the residual (for 
small scale spatial variation). Figure [5] suggests clear large scale spatial structure depending on 
latitude. Thus, we modeled the mean structure through simple harmonic regression: 


m(L) = oo + Oi cos(L x 7r/90°) + 02 sin(L x 7r/90°). (8) 

We considered two cases: one is with a constant (unknown) mean (that is, Oi = 02 = 0 in ([8])), 
and the other is given by ([8]). In both cases, we hrst estimate the mean structure using regression 
and then work with the residual to £t the covariance structure. 


4-2. Example I: horizontal directional sampling design for prediction 

Figure [6](a) shows the map of residual after subtracting the constant mean and Figure [6](b) 
shows the empirical semi variogram of residuals. The semi variogram clearly shows negative co- 
variances for large distances. To save computational burden, we randomly selected 600 locations 
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near the red region for parameter estimation, and selected 200 locations over the blue region for 
prediction as in Figure |6]^a). We repeated this process 100 times and for each time, all locations 
are randomly sampled and thus different. We compared the covariance models, MC, MG, and 
the convex sum model, C, with G'^-Wendland covariance function ([5]) and cosine function. As 
shown in Table [21 MG gives poor model fit and prediction. Although MG has the largest max¬ 
imum log-likelihood value, it leads to poor prediction compared to G. This is expected as the 
Matern model is not able to produce negative correlations. Similarly to the simulation example 
of Section 13.21 G is the best model in terms of prediction. 

Figures a) and (b) show boxplots of differences of AE and GRPS values between MG and 
G, displayed against minimum great circle distance between a prediction location and its nearest 
sampling location. For all distances, G outperforms MG. Moreover, the differences of two pre¬ 
diction errors between two models increase as minimum distances between prediction locations 
and their nearest sampling locations increase. On the other hand. Figure IHl^a) presents residual 
helds after removing mean structure by using simple harmonic regression depending on latitude 
as in ([8]). We htted (PAWendland covariance functions with the Euclidean distance (WG) and 
the great circle distance (WG) in addition to the covariance models considered previously with 
constant mean structure. For WG, the covariance function is dehned by 

= a^{l + {tT)/{cR) + 

where r > 6, c > 0, and t > 0. Sampling and prediction locations remained the same as 
in the previous example. From Table [3l all models have comparable maximum log-likelihood 
values except MG. This may be due to the fact that the residual field seems smooth. Regarding 
prediction, G and WG have better performances than MG and WG, respectively. Note that WG 
has much larger sample standard deviations of estimates for support and shape parameters than 
WG. 

When comparing the two mean structures, prediction errors for the case of a constant mean 
are larger than those for the case of the mean given by (|H]). Since the geopotential height data 
mainly show large scale, smooth, variation depending on latitude, the mean structure using 
simple harmonic regression resulted in improved prediction. Overall, the covariance functions of 
the class performed better than those of the class ‘Fd+i regardless of mean structures. 
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4-3. Example II: vertical directional sampling design for prediction 

We entertain the same set of covariance models with the mean strnctnres as Section 14.21 
However, we changed a sampling design for prediction. We randomly selected 600 locations 
where longitnde is less than 0 ° for parameter estimation, and selected 200 locations where that 
is greater than 0° for prediction as in Fignre |9]^a). The empirical semi variogram in Fignre[9](b) 
shows that there is not mnch non-negative covariance valnes, nnlike the previons example. 

With a constant mean strnctnre. Table H] shows similar results as in Section 14.21 All models 
except MG are similar in terms of maximum log-likelihood values. Both C and WG have better 
performance in terms of prediction than MG and WG, respectively. Although results in Table 0] 
show that models G and WG do not outperform MG and WG as signihcantly as in Section 
there still exist some improvements in terms of prediction with the models dehned with the great 
circle distance. From Figures [TUf al and (b), we observe that WG has smaller AE and GRPS 
values than WG for prediction locations who are relatively far away from their nearest sampling 
locations. 

When we consider the mean structure in (IH]), there is no signihcant difference in terms of 
prediction errors between G^-Wendland models using the great circle distance and the Euclidean 
distance. However, convex sum and G^-Wendland models using the great circle distance perform 
better than Matern covariance model using the Euclidean distance in prediction. It is expected 
that the vertical directional design has smaller prediction errors than the horizontal directional 
design in Section 14.21 (see Tables [2] and |4]) due to large scale variation depending strongly on 
latitude. 

5. Discussion 

We have considered several classes of isotropic covariance functions with either the great circle 
distance or the Euclidean distance and compared them in terms of parameter estimation and 
spatial prediction. We have shown that when the true spatial range is large, the prediction 
performance of covariance models dehned on the sphere using the great circle distance (that is, 
'4>{6) on Sf) is better than the functions projected from the Euclidean space in spatial prediction. 
Moreover, when data showed signihcantly negative correlations at large distance lags, isotropic 
covariance models in the class $3 are not adequate and there is substantial diherence between 
covariance models from the classes 4/2 and $3 in prediction. In geopotential height data set, we 
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showed that the distortion of the Euclidean distance may lead to poor prediction for prediction 
locations that are relatively far away from sampling locations. 
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Table 1: (Simulation example on 5^) Sample means and standard deviations of parameter es¬ 
timates, maximum log-likelihood values, and prediction errors for each model (100 cases). For 
model C, A G (0,1) is a weighting parameter. 


Model 

MC 

MG 

C 

H 


6.654(1.237) 

3.260(0.448) 

26.033(7.163) 

7.496(1.711) 

d or A 

7918.459(1669.169) 

49423.925(7016.223) 

0.999(0.006) 

58.499(56.164) 

u or 13 

1.037(0.054) 

0.500(0.000) 

1.820(0.038) 


Max.loglik 

356.505(11.217) 

298.233(6.923) 

354.384(10.533) 

-430.486(38.393) 

RMSE 

1.123(0.179) 

0.975(0.085) 

0.874(0.126) 

1.961(0.369) 

MAE 

0.903(0.177) 

0.749(0.065) 

0.717(0.118) 

1.751(0.337) 

CRPS 

0.639(0.096) 

0.573(0.058) 

0.514(0.055) 

1.166(0.213) 


Table 2: (Data example I - constant mean) Sample means and standard deviations of parameter 
estimates, maximum log-likelihood values, and prediction errors for each model (100 cases). For 
model C, c G (0, tt] is a support parameter and r > 6 is a shape parameter. 


Model 

MC 

MG 

C 


1.410(0.087) 

0.800(0.115) 

2.240(0.377) 

d or A 

2128.225(114.514) 

136377.194(18394.328) 

0.811(0.108) 

z) or c 

2.605(0.072) 

0.500(0.000) 

2.931(0.084) 

f 



7.826(0.302) 

Max.loglik 

2167.401(20.798) 

1413.049(11.086) 

2160.817(20.348) 

RMSE 

3.972(0.137) 

6.939(0.086) 

3.292(0.186) 

MAE 

3.853(0.142) 

6.840(0.090) 

3.168(0.193) 

CRPS 

3.427(0.145) 

6.767(0.089) 

2.650(0.196) 
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Table 3: (Data example I - simple harmonic regression depending on latitude) Sample means 
and standard deviations of parameter estimates, maximum log-likelihood values, and prediction 
errors for each model (100 cases). For models C and WG, c e (0, vr] is a support parameter and 
r > 6 is a shape parameter. For model WG, c > 0 is a support parameter. 


Model 

MC 

MG 

C 

WG 

WC 

<72 

0.776(0.076) 

0.100(0.001) 

0.820(0.053) 

0.706(0.036) 

0.770(0.139) 

d or A 

1770.346(121.073) 

10754.167(91.314) 

0.798(0.056) 



P or c 

2.699(0.088) 

0.500(0.000) 

2.887(0.074) 

2.873(0.020) 

9.794(22.488) 

f 



9.426(0.225) 

9.237(0.198) 

29.986(68.228) 

Max.loglik 

2171.231(21.331) 

1242.183(6.178) 

2168.634(6.178) 

2168.689(20.324) 

2168.653(20.271) 

RMSE 

2.441(0.137) 

2.351(0.057) 

2.177(0.085) 

1.988(0.074) 

2.184(0.357) 

MAE 

2.314(0.142) 

2.217(0.065) 

2.047(0.090) 

1.856(0.079) 

2.053(0.365) 

CRPS 

1.952(0.135) 

2.100(0.064) 

1.677(0.084) 

1.493(0.072) 

1.689(0.351) 


Table 4: (Data example II - constant mean) Sample means and standard deviations of parameter 
estimates, maximum log-likelihood values, and prediction errors for each model (100 cases). For 
models G and WG, c G (0, tt] is a support parameter and r > 6 is a shape parameter. For model 
WG, c > 0 is a support parameter. 


Model 

MC 

MG 

C 

WG 

WC 

<72 

2.062(0.278) 

0.311(0.005) 

3.621(0.315) 

3.545(0.302) 

3.249(0.208) 

d or A 

1220.475(158.687) 

10527.394(87.840) 

0.934(0.084) 



U OT C 

3.412(0.229) 

0.500(0.000) 

2.939(0.168) 

3.026(0.061) 

3.018(0.957) 

f 



8.539(0.648) 

8.688(0.377) 

8.783(2.826) 

Max.loglik 

1959.170(40.787) 

921.716(12.424) 

1923.248(31.613) 

1923.297(31.593) 

1924.991(32.035) 

RMSE 

1.167(0.079) 

1.610(0.063) 

1.107(0.118) 

1.042(0.090) 

1.068(0.097) 

MAE 

0.929(0.075) 

1.340(0.063) 

0.855(0.105) 

0.804(0.086) 

0.832(0.090) 

GRPS 

0.606(0.054) 

1.167(0.059) 

0.558(0.064) 

0.527(0.053) 

0.541(0.057) 
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Figure 1: (Simulation example on iSj^) A realization of sampling locations (o) and prediction 
locations (x). 



0 (radian) 


9 (radian) 


0 (radian) 


0 (radian) 


Figure 2: (Simulation example on iSj^) MAE (a) and mean CRPS (b) averaged over 100 replica¬ 
tions from the two models displayed against prediction locations (0) when the true spatial range 
a = 27r. (c) and (d) are the same as (a) and (b), except that a = vr/d. Triangles and circles 
represent the values of prediction errors for the exponential models using great circle and chordal 
distances, respectively. 
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Figure 3: (Simulation example on iS|.) (a) A realization of residual fields. For (a), sampling 
locations (o) and prediction locations (x). (b) Empirical semi variogram values for selected 
locations versus the great circle distance. For (b), dotted line represents sample variance. 
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Figure 4: (Simulation example on iS|.) Boxplots of differences of AE (a) and CRPS (c) values 
from MC and MG, displayed against minimum great circle distance between a prediction location 
and its nearest sampling location, (b) and (d) The same as (a) and (c), except that selected 
models are MG and G. For (a)-(d), red circles represent average values in each bin. 
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Figure 5: (Data example I) Square root of the geopotential height at level 500 hPa 
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Figure 6: (Data example I) (a) A realization of residual fields after subtracting the constant mean. 
For (a), sampling locations (o) and prediction locations (x). (b) Empirical semi variogram values 
for selected locations versus the great circle distance. For (b), dotted line represents sample 
variance. 
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Figure 7: (Data example I - the constant mean) Boxplots of differences of AE (a) and CRPS (b) 
values from MC and C, displayed against minimum great circle distance between a prediction 
location and its nearest sampling location. Red circles represent average values in each bin. 
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Figure 8: (Data example I) (a) A realization of residual helds after removing mean structure 
through simple harmonic regression depending on latitude. For (a), sampling locations (o) and 
prediction locations (x). (b) Empirical semi variogram values for selected locations displayed 
against the great circle distance. For (b), dotted line represents sample variance. 
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Figure 9: (Data example II) (a) A realization of residual fields after subtracting the constant 
mean. For (a), sampling locations (o) and prediction locations (x). (b) Empirical semi variogram 
values for selected locations displayed against the great circle distance. For (b), dotted line 
represents sample variance. 
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Figure 10: (Data example II - the constant mean) Boxplots of differences of AE (a) and CRPS (b) 
values from WG and WC, displayed against minimum great circle distance between a prediction 
location and its nearest sampling location. Red circles represent average values in each bin. 
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